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Abstract 

On January 14, 2005, ESA’s Huygens probe separated from NASA’s Cassini 
spacecraft, entered the Titan atmosphere and landed on its surface. As part of 
NASA Engineering Safety Center Independent Technical Assessment of the 
Huygens entry, descent, and landing, and an agreement with ESA, NASA 
provided results of all EDL analyses and associated findings to the Huygens 
project team prior to probe entry. In return, NASA was provided the flight data 
from the probe so that trajectory reconstruction could be done and simulation 
models assessed. Trajectory reconstruction of the Huygens entry probe at Titan 
was accomplished using two independent approaches: a traditional method and 
a POST2-based method. Results from both approaches are discussed in this 
paper. 


INTRODUCTION 

On January 14, 2005, the European Space Agency’s (ESA) Huygens probe 
separated from NASA’s Cassini spacecraft, entered the Titan atmosphere and landed on 
its surface. 1 "’ This probe used a multiple parachute system to enable atmospheric 
measurements to be recorded during the probe’s over two hour descent to the surface. 
Digital images, radar altimetry, accelerometer data, and Earth-based radio telescope 
observations were also gathered during the entry, descent, and landing (EDL). 4 " 9 Figure 1 
illustrates the Huygens probe’s EDL profile. After atmospheric interface at 1270 km 
above the surface, the probe decelerates to around Mach 1.5 at pilot parachute deploy. 
This deploy event (designated TO) is triggered by a sequence of time and acceleration 
conditions and is the epoch time for all subsequent events and most data sets generated 
during the parachute descent phase. Three parachutes were used in the Huygens probe 
system: (1) a 2.5-sec Pilot parachute; (2) a 15-minute Main parachute; and (3) a 2.5-hour 
Drogue parachute. Data taken during the descent was relayed through the Cassini 
spacecraft as it flew by Titan. 
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NASA Langley was involved in the pre-entry EDL analyses of the Huygens 
probe. Analyses were conducted under the auspices of the NASA Engineering Safety 
Center’s (NESC) Independent Technical Assessment (ITA) of the Cassini/Huygens probe 
EDL at Titan. 10 A Program to Optimize Simulated Trajectories II 11 (POST2)-based 
trajectory simulation was developed that included models of the probe aerodynamics, 
parachute trigger logic and drag models for the Pilot, Main, and Drogue parachutes. As 
part of the agreement with ESA, NASA provided results of all analyses and presented 
findings to both the Cassini and Huygens project teams. In return, NASA was provided 
the flight data from the probe so that trajectory reconstruction could be done and 
simulation models assessed. NASA has completed similar assessments of flight data to 
improve simulation models (such as capsule aerodynamics, parachute aerodynamics, and 
atmospheric density and winds) for the Mars Pathfinder EDL, Mars Odyssey 
aerobraking, and Mars Exploration Rovers EDL. 12 14 ESA had a team that provided the 
official project reconstruction of the Huygens EDL trajectory (the Huygens Descent 
Trajectory Working Group, or DTWG). 3 The main objective of this NESC sponsored 
activity was to reconstruct the Huygens Probe data to improve NASA’s aerodynamics, 
atmospheric density and winds, and parachute performance models. The NASA NESC 
ITA team also analyzed the flight data to provide an independent trajectory 
reconstruction. The results of these analyses were provided to the DTWG and the 
Huygens Data Analysis Working Group (DAWG), with interaction between NASA, 
DTWG, and DAWG to discuss any differences in the reconstructed trajectory. 
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BACKGROUND 

Trajectory reconstruction of the Huygens entry probe at Titan was accomplished 
using two independent approaches. One was a traditional approach to trajectory 
reconstruction that integrates the accelerometer measurements directly during the entry 
phase (forward) and integrates the hydrostatic equation using the temperature and 
pressure measurements to determine altitude and velocity (backwards). The other 
approach used a Kalman filter module developed for reconstruction in conjunction with 
the POST2-based simulation. The POST2-based reconstruction uses accelerometer 
measurements to adjust an estimated state using a POST2-based simulation developed 
prior to entry to support EDL analyses and design. Results from both approaches are 
discussed in this paper. 

The ESA Huygens probe science teams provided several sets of data. These data 
included acceleration measurements throughout EDL, pressure and temperature after 
heatshield jettison, and radar altimetry for the last 40 kilometers. Additionally, radio 
telescopes at Earth received the Huygens signal throughout the descent and provided an 
estimate of wind velocity. The Descent Imager (DI) team also contributed an estimate of 
vertical velocity based on sequentially captured images. 

The main emphasis of the Huygens probe reconstruction is to evaluate the 
simulation models used prior to probe entry at Titan with actual flight data. Another 
objective of the Huygens reconstruction was to compare the NASA derived trajectory and 
the Huygens DTWG solution to the trajectory profile. 


RECONSTRUCTION APPRO ACHS 

Traditional Method 

Using the traditional method, the trajectory reconstruction of the Huygens probe 
trajectory parameters from entry interface to the surface of Titian is separated into two 
parts. First, the probe trajectory on the Drogue parachute (lasting about 8650 s) is 
addressed. Next the probe entry phase trajectory (the main deceleration pulse), leading to 
the parachute phase and culminating just prior to the Drogue deploy (lasting about 350 s) 
is reconstructed. 

The reconstruction of the trajectory during the Drogue phase of the Huygens 
probe descent to the surface is obtained using the pressure, temperature, and acceleration 
data. In essence, the method involves integrating a form of the hydrostatic equation to 
obtain altitude, and applying the aerodynamic force equation to get total relative velocity 
magnitude. This approach is taken to circumvent the problems associated with integrating 
the accelerometer data during long time intervals at near terminal velocity conditions. 

The second half of the reconstruction was done in a more classical sense. It 
involved the hypersonic portion of the flight starting from entry interface. By using the 
acceleration data along with the navigation team provided state at entry interface the 
trajectory was reconstructed by integrating the accelerations. After the trajectory was 
reconstructed the atmosphere could be backed out using the aerodynamic database. 
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Atmosphere Reconstruction 

The accelerometer data from an entry probe during descent is the principle data 
used to obtain the free-stream density. The aerodynamic force equation along the x-body 
axis can be expressed as, 


P = 


2ma x 

v;c A s 


where the probe mass, m, reference area, S, and the axial aerodynamic force coefficient, 
Ca are known a priori, and v r is the probe velocity relative to the atmosphere. Given the 
reconstruction trajectory state variables (i.e. probe position, velocity, and orientation as a 
function of time); density mapping is relatively straight forward, except for the Ca, which 
is also a function of density (along with other parameters). This transcendental nature of 
the density equation above can be overcome by a simple iterative technique. 

The aerodynamic coefficient for the hypersonic phase, Ca, is determined using the 
preflight database created by NASA Langley. This database requires the current values 
for the center-of-gravity location, Knudsen number, Mach number, relative velocity; 
angle of attack and sideslip angle are assumed zero. During the parachute/drogue phase, 
the composite CaS coefficient is determined as a linear combination of probe and 
parachute drag. That is, for the total CaS term the product of the approximate drag 
coefficient, Cd (.8) and the heat shield reference area (5.73), or the descent module 
reference area (1.29) was added to the parachute drag (a function of Mach Number and 
Reynolds Number) times its reference area. 

Drogue-Phase Equations of Motion (Hydrostatic equation) 

The equation for the change in pressure due to a change in altitude of a fluid, such 
as an atmosphere, is known as the hydrostatic equation, and is given as: 


dp 

dh 


= -g P 


( 1 ) 


Using the chain rule, substituting in equation (1) and solving for the altitude time 
derivative gives: 


dh p 
dt -g p 


( 2 ) 


Putting the density (from the equation of state using the measured temperature and 
pressure) and gravity into equation 2, gives a transcendental equation that is a function of 
the measured quantities that can be evaluated and integrated to get the vertical velocity 
component and the altitude time history. 

dh R(R 0 + h) 2 T / p ) 
dt -^C f W MM [pj 


( 3 ) 
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This equation is integrated using the Ordinary Differential Equation (ODE) solver 
that is built into MatLab. It requires that the pressure, temperature, mean molecular 
weight, and pressure derivative be given as a function of time. The gravity and 
compression factor are functions of altitude but since an ODE solver is being used, which 
uses the previous time point to calculate the derivative to get to the next point, the 
altitude is known. Since the time of surface impact is known, the integration is done from 
the surface up. In summary, from the measurements of pressure and temperature as a 
function of time, two of the six trajectory state variables (h, h) as a function of time are 
obtained. 

POST2-based Method 

A six degree-of-freedom atmospheric entry and three degree-of-freedom 
parachute descent trajectory of the Huygens probe was simulated in Program to Optimize 
Simulated Trajectories II (POST2). 11 POST2 is a generalized point mass, discrete 
parameter targeting and optimization trajectory program. POST2 has the ability to 
simulate three and six degree-of-freedom (3DOF and 6DOF) trajectories for multiple 
vehicles in various flight regimes. POST2 also has the capability to include different 
atmosphere, aerodynamics, gravity, propulsion, parachute and navigation system models. 
Many of these models have been used to simulate the entry trajectories for previous 
missions (i.e., MER, Genesis, Mars Pathfinder) 15 " 17 as well as current and planned NASA 
missions (Stardust, Mars Phoenix, and Mars Science Laboratory). 18 ' 19 A variety of 
system studies have been conducted and their atmospheric trajectories simulated in 
POST2 including aerocapture at Titan, Neptune and Venus. 20 

Exploiting the modular nature of the POST2 program by adding mission specific 
models in concert with the existing POST2 architecture allows for the development of 
higher fidelity, mission specific simulations. These simulations usually have integrated 
mission-specific engineering models and flight software that have been tested and 
validated prior to the actual mission. These simulations support design, development, 
testing, and operations of vehicles for particular missions. These POST2-based mission 
specific simulations have been used operationally for day-of-entry flight software 
parameter determination and prediction of mission metrics (such as touchdown footprint) 
as well as aerobraking orbit prediction and assessment. Simulation complexity varies 
from first-order trades (e.g. parachute size and deployment conditions, terminal descent 
engine size, etc.) to all-up Monte-Carlo simulations for flight operations. 

POST2 was used to simulate the Huygens entry trajectory into Titan. 10 The 
POST2-based flight simulation incorporated several models specific to the Huygens 
probe entry: a 6DOF aerodynamics model; Titan’s gravity and Titan-GRAM atmosphere 
(including wind) models; attitude inputs and initial states; trigger criteria, inflation, and 
drag models for the Pilot, Main, and Drogue parachutes; as well as vehicle geometric 
parameters. Note that version 1.0 of the Titan-GRAM atmospheric model was 
implemented into POST2, with updates from Cassini measurements of Titan (the TO and 
TA atmospheric profiles). This simulation was used to produce single trajectory data and 
was an integral element of the Monte Carlo analyses discussed in Ref. 10. 
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Extended Kalman Filter P0ST2 Module 

An extended Kalman filter (EKF) module has been developed for and included 
into POST2. This module was developed and integrated into the POST2 software to 
facilitate trajectory reconstruction using POST2-based mission specific simulations. The 
general nature of POST2 inputs was retained for this module. That is, the state to be 
estimated can be input as any POST2 input and the observation can be any POST2 
outputs. While the general POST2 software architecture was retained, separate files of 
observations and their associated weightings versus simulation time are required for use 
with this EKF module. The utility of integrating this EKF function into POST2 is to 
allow more rapid setup and execution of trajectory reconstruction runs using the same 
simulation that has been tested and validated for that particular mission. 

The theory and equations defining this module are well described in the 
literature. 21 ' 22 A summary of the method implemented for the POST2 module is as 
follows: 

1) Define initial covariance and estimated state (Po and Xo) 

2) Read observations (O) and their weights (W = R 1 ) from files 

3) Integrate to the next observation time: 

X = F(X,t) and F = A* P + P * A T + Q- K* R* K T 

where P is the propagated state covariance matrix; A is the matrix of state derivative 

dX 

with respect to time (A) sensitivity to the state ( — ); Q is the system noise covariance 

dX 

matrix (where the random noise is chosen to apply to every state); matrix K is from the 
last update; and R (the measurement noise covariance) is set to the inverse of the 
observation weightings matrix (W). 

4) At the observation time, calculate: 

Y = O - G H = — K = P*H T *(H*P*H T -R)~ l 

dX 

where O is the observation matrix (from file), G is the matrix of calculated observation 
values from the simulation and K is the Kalman gain. 

5) Next update the state: 

X = X + K*Y 

6) And update the covariance: 

P = (I-K*H)*P*(I-K*H) T + k*r*k t 

where I is the identity matrix, X is the estimated state, and P is the updated covariance. 
This process is continued with POST2 integrating the state and covariance between 
observation updates. 
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Huygens EDL Measured Datasets 

The Huygens entry, descent, and landing datasets described below were provided 
to NASA as part of an agreement with ESA for analyses NASA performed prior to entry. 
These datasets include measurements of atmospheric pressure, temperature, composition, 
and zonal winds, vehicle acceleration, and vehicle altitude. Various Principal 
Investigators were responsible for the instruments that provided this data; the various 
groups responsible for this data are identified, along with details about their instruments, 
in Refs. 1 and 2. 

Pressure 

The time history of the reduced pressure was obtained from ESA as described 
above. Figure 2 (top panel) shows the free-stream pressure from the onboard pressure 
transducer as a function of time where the reference time is entry interface. (During the 
day of entry, Jan. 14, 2005, t E i = 09:05:52.523) 




Figure 2: Upper Panel: Free-stream pressure provided by ESA. 

Lower Panel: Pressure derivative using a 51 -point sliding window with a third 
degree polynomial evaluated at the mid-point. 


Derivative of Pressure 

The time derivative of pressure is needed for the calculation of altitude, see 
equation (3). The method used to obtain the pressure derivative includes fitting a third 
degree polynomial (using least-squares) to a 51 -point sliding window of pressure data. 
That is, at any time point, 25 pressure data to the left and right are fitted to a 3 rd -degree 
polynomial from which a derivative in obtained for the time point. One point to the right 
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is added, while one point to the left is dropped from the original data set and the process 
is repeated. This method clearly leaves 25 points at the beginning and at the end of the 
data set without a derivative. To handle these end points, the window size is reduced to 5 
points, leaving only two points without a derivative. Figure 2 (lower panel) shows the 
results of the calculations for the time derivative of pressure, dp/dt. The derivative is 
relatively smooth except at the end data points, as expected. The left end of the graph has 
a little extra noise providing a few points with an unreliable derivative. On the ground, 
the derivative goes to zero and this “corner” impacts the last few derivatives, as seen in 
the graph. Neither end produces difficulty in the integration process since the points can 
be readily omitted. 

Density 

The density is obtained from the calibrated measurements of pressure and 
temperature using the equation of state. In addition, the mean molecular weight and the 
compression factor (also taken from data provided through ESA, but not shown here) are 
used in the calculation of density. Figure 3 shows the density results for both with and 
without (Cf=l) compression factor. As expected, the difference is less than 3.5% and the 
larger differences occur at lower altitudes, corresponding to larger values of time, as 
shown in Fig. 3. 



Figure 3: Density time history with and without compression effects. 


Altimeter 

Two separate independent altimeters, labeled “A” and “B”, were onboard the 
Huygens probe. Figure 4 shows the data for both of these instruments after the Principal 
Investigator has applied calibration factors. The data produced by both altimeter 
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instruments lay on top of one another, although the altimeter B data starts at a higher 
altitude. 




Figure 4: Top Panel : Altimeter measurements from A and B instruments. 

Lower Panel'. Derivative of altimeter data using a 101 -point moving window. 

Accelerometer 

The axial acceleration data set from the Huygens Servo Accelerometer, used in 
the reconstruction is shown in Fig. 5. Note that Piezo Accelerometer data in all three axes 
were received, but evaluation of this dataset deemed it unusable (the axial component 
only reached 100 m/s 2 , well below the values expected and returned from two 
independent sources, and forcing the lateral measurements to be rejected as well). 
Examination of the Servo axial accelerometer data before entry interface indicated that 
there is a bias in the data of about 22.8 p-g in the data. In Fig. 5 the top panel shows the 
data from entry interface to near probe touchdown on the surface of Titan. The lower 
panel shows an expanded view of the data. In this panel the results of a smoothing 
process is shown along with the original data. Smoothing was done on data beyond 330 s, 
or right after parachute inflation. The smoothing process consists of a median process for 
a sliding window, similar to what was done for the calculation of dp/dt and dh/dt, 
discussed earlier, except a median is obtained instead of a curve fit. That is, a window 
size of 201 data points is selected and a median of the data set is determined. The time 
and median value are recorded at the middle of the data and subsequently, a data point is 
added to the right and one deleted to the left and the process is repeated. Again, because 
of the moving window the first and last 100 points of data in the original data set remain 
not smoothed. The first 100 points are used to merge with the original data set, while the 
last 100 points are not smoothed, as can be seen in Fig. 5. The spike in the acceleration 
data at approximately 1100 seconds corresponds to drogue deploy which represents real 
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data and as such the peak must be preserved so as to not lose the information. Thus, in 
this time interval smoothing was not attempted. Similarly, the hypersonic portion of the 
acceleration time history before main parachute deploy, i.e. before approximately 330 
seconds, is also not smoothed so as not to lose the peak acceleration data. In summary, 
the blue line in Fig. 5 represents the acceleration data used in the calculations reported in 
the next sections. 
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Figure 5: Upper Panel: Huygens axial accel from Entry Interface to Titan’s surface. 
Lower Panel: Expanded view of axial acceleration showing smooth and actual accel data. 


Expanded View 



RESULTS 

Traditional Method 

Figures 6, 7, and 8 show the atmosphere results for the hypersonic reentry phase 
of the mission (denoted “Hyper” in the legend), namely density, pressure, and 
temperature as a function of altitude, respectively. Included on each graph is the Huygens 
model atmosphere generated by Yelle and the Titan-GRAM Atmosphere used in 
simulations. Also included in the figures is the lower atmosphere reconstructed density 
using the equation of state (denoted “Hydro” in the legend) as well as the measured 
pressures and temperatures. 
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Figure 6: Reconstructed atmospheric density as a function of altitude. 



Figure 7: Reconstructed atmospheric pressure as a function of altitude. 
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Figure 8: Reconstructed atmospheric temperature as a function of altitude. 

The atmosphere state (density, pressure, and temperature) graphs also include the 
lower altitude measurement data discussed earlier for completeness. At the lower altitude 
(approximately less than 500 km) there is an excellent match between the Yelle model 
and the measurements, and an even better match with the Titan-GRAM. At altitudes 
larger than about 500 km, both the density and pressure give higher values at the 
exosphere base resulting in temperature that agrees with the model. Between about 500 
and 800 km, the model appears to diverge slightly, resulting in temperature differences of 
only about 30 K. For the entire atmosphere profile the Titan-GRAM model agrees more 
closely to the reconstruction than the Yelle model. To facilitate the comparison of the 
reconstructed density to the models a comparison chart showing the percent difference 
between the reconstruction and the model was created and is presented as Fig. 9. 



Figure 9: Percent Difference Between Reconstructed Density and Model Densities 
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Drogue Trajectory Reconstruction 

Figure 10 shows the altitude results from the integration of equation (3) starting at the 
ground and ending near the time when the heat shield is jettisoned. The change in the 
slope of altitude after about 1100 s shown in the graph is the drogue deployment. 
Included on the graph are both of the altimeter measurements taken at altitudes starting at 
about 40 km. 



Figure 10: Reconstructed altitude compared to both altimeter measurements. 

Figure 11 shows a more detailed comparison between the reconstructed altitude and the 
altimeter-A measurements (as mentioned above, both altimeters-A and B provided nearly 
the same data, so only altimeter-A will be shown). In addition, the corresponding 
residuals are plotted to show the detail differences. Clearly, as seen in the lower panel for 
Fig. 11, there is a systematic difference between the reconstructed altitude and the 
altimeter measurements for this altimeter (the same differences of about the same 
magnitude apply to altimeter-B). 
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Figure 1 1 : Upper Panel : Altimeter A comparison to the reconstructed altitude. 

Lower Panel'. Altitude measurement residuals with the RMS. 

Descent Phase Comparisons 

Figure 12 shows a comparison of the reconstructed altitudes obtained by the 
DTWG with the altitudes developed using other reconstruction methods. There is 
excellent agreement with the DTWG reconstructed altitudes using the traditional 
reconstruction method described in this report (the POST2-based method results are 
discussed in following sections). The two curves lie almost on top of each other (about 
300 m difference at the upper altitudes to about 10 m near the surface). This result is 
mainly due to the similarity in reconstructed methods. Namely, the DTWG method uses 
an approximate integral, while numerical integration is used herein. Clearly, the DTWG 
method does not match the altimetry measurements. 
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Figure 12: Reconstructed altitude comparison with RADAR 

Figure 12 also shows a similar comparison of the reconstructed altitudes obtained 
by Carlo Bettanini, a member of the Huygens Atmospheric Structure Instrument team 
(personal communication, 20 September 2006). There is good agreement between these 
results also. Once again the reconstruction does not match the radar altimetry data; note 
that Bettanini’s method does not use the altimetry measurements to reconstruct altitude. 
His method relies upon integration of the pressure measurements combined with 
accelerometer data from the ground upward. There are differences at higher altitudes 
(above 105 km) where Bettanini altitudes are lower in value than what is generated using 
the methods in this report. In summary, the altitude results using the traditional method 
match the DTWG and Bettanini reconstructed altitudes. However, when compared to the 
radar altimetry data, a systematic error is readily evident. One possible solution is to 
introduce a scale factor on the radar altimetry data. 

Altitude Error Bound Estimation 

A useful metric for understanding the altitude variation is to estimate the error in 
the altitude reconstruction, i.e. determine error bars or bounds. Since the altitude is found 
by integrating the measured data it stands to reason that a main source of uncertainty in 
the altitude would come from the uncertainty in the measurements. The measurements of 
the pressure, temperature, and mean molecular weight each have their own associated 
uncertainty. To find the extent of the error bounds for velocity the two extreme cases of 
measurement uncertainty were taken, that is all were either high or low. That is, changing 
each measurement by adding its uncertainty to it for the high case and subtracting for the 
low case. After adjusting the measurement, the integration was done again to find the 
upper and lower bounds of the error bar on altitude. Because of the large scale of the 
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altitude, 0 to 150 km, the error bars are difficult to visualize. Therefore it is easier to 
show the size of the error bounds as a function of altitude rather than the actual bounds 
on either side of the altitude history. This error is shown in Fig. 13, where it can be seen 
that the error bounds are relatively small in comparison to the actual altitude. 



Figure 13: Altitude Error Bound Size as a function of Altitude 


POST2-bcised Method 

The trajectory reconstruction runs using POST2 focused on estimating entry 
aerodynamics, parachute phase parameters, and the vehicle state (position and velocity) 
throughout the entire EDL trajectory. Although the simulation results using the JPL- 
determined last estimated nominal entry state generally follows the flight data, a Monte 
Carlo analysis was used to try to determine the initial state (inertial position and velocity 
vectors) that best fits the flight data through the entry phase (where the main deceleration 
pulse occurs). The last JPL delivered initial state covariance was used to generate 10000 
dispersed states. Each case was simulated and the root sum square (RSS) was taken of the 
acceleration observation residual as the simulation progressed. The initial state from the 
smallest value cases was used for subsequent reconstruction efforts. Table 1 compares the 
nominal best estimated entry state generated by JPL (along with the 3-sigma bounds) and 
the initial state determined by the Monte Carlo process (and the differences between 
them). Note that all times in the following sections are from entry interface and not Pilot 
parachute deployment (or TO). 
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Table 1. Comparison of JPL Last Initial State Estimate and Monte Carlo Determined 
Initial State 



JPF last 
estimate 

3-o bounds 

Monte Carlo 
determined 

Difference 

Inertial X Position (km) 

-3785.053 

94.718 

-3810.021 

24.968 

Inertial Y Position (km) 

366.623 

24.508 

367.226 

-0.603 

Inertial Z Position (km) 

-568.429 

10.214 

-566.418 

-2.011 

Inertial X Velocity (km/s) 

5.704492 

0.009794 

5.701933 

0.002559 

Inertial X Velocity (km/s) 

1.918924 

0.001877 

1.918745 

0.000179 

Inertial X Velocity (km/s) 

0.390306 

0.001727 

0.389695 

0.000611 


Entry Aerodynamics 

The focus of this section is to validate the entry aerodynamics model from the 
pre-entry simulation using the flight data. For this portion of the analysis, the Titan- 
GRAM atmosphere model is assumed correct and all of the error in the measured 
accelerations is due to errors in the aerodynamics model. The POST2-based Huygens 
EDL simulation was setup to include aerodynamic uncertainties for Monte Carlo 
simulations done prior to probe arrival at Titan. The aerodynamic database includes the 
uncertainty amount for the free molecular, hypersonic and supersonic flight regimes. For 
the reconstruction analyses, the uncertainty in axial force coefficient was estimated since 
only the axial accelerations were measured. Note that these were entry phase runs only, 
hence no reconstructed trajectory data is generated beyond the Pilot parachute mortar 
firing event (aka TO event) at about 270 sec. 

As noted previously, the Kalman filter implementation includes system noise 
covariance input. Through the state covariance, this noise covariance can be used to 
increase the reconstructed state response to the measurements. The peak deceleration 
region during entry is shown for the reconstruction case with the noise covariance set for 
a rapid filter response to the measurements in Fig. 14. This “Fow Noise” (or FN) setting 
allows an assessment of the entry aerodynamic model used prior to Huygens probe entry 
at Titan by assuming that all of the error during the entry is from aerodynamic 
misprediction; whereas a “High Noise” (or HN) setting does not allow the filter to change 
the estimate as rapidly. As seen in this figure, the acceleration for the FN case matched 
acceleration points closely as desired. The maximum error during this peak deceleration 
period was +/- 0.6 m/s 2 (or less than 0.5%); most of the error in this region was less than 
0.2 m/s 2 (or 0.2%). 
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Figure 14. Comparison of Flight Data to Entry Aerodynamics Estimate Reconstruction 
Run - Peak Deceleration Region - Low Noise 

These reconstruction runs were estimating the adjustment from the nominal 
aerodynamics model necessary to meet the flight data. The actual states estimated were 
the aerodynamic dispersion parameters for Ca in each flight regime (free molecular, 
hypersonic and supersonic). These parameters varied between +/- 1 and applied the 
aerodynamicist-deflned uncertainty to the coefficient in each regime. For example, in the 
free molecular flight regime an uncertainty parameter of +1 corresponds to a multiplier of 
5% on Ca; that is, the nominal Ca was multiplied by 1.05 for an uncertainty parameter of 
+1 in the free molecular flight regime. 



Figure 15 Estimated Entry Axial Force Aerodynamic Uncertainty Parameters -Low Noise 
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The estimated uncertainty parameters for the LN reconstruction case are shown in 
Fig. 15. Note that, the aerodynamics model allowed for overlap between the flight 
regimes, hence there are regions of the plot where two parameters are active at the same 
time. Otherwise, the unused uncertainty parameter is a constant value. The key result 
from this figure is the fact that even with all the error during entry attributed to the 
aerodynamics the estimated solution never exceeds the 3 -sigma uncertainty bounds 
(never greater than +1 or less than -1) established by the aerodynamicists prior to entry. 
This significant result indicates that the entry aerodynamics model captured the 
aerodynamics of the Huygens entry capsule throughout all three flight regimes. 



Figure 16 Estimated Entry Aerodynamic Axial Force Coefficient Values 


Figure 16 shows the actual axial coefficient values corresponding to the estimated 
uncertainty parameters throughout entry for both the HN and LN versions of this 
reconstruction case as well as the simulation only case as a function of relative velocity. 
The 3-0 axial force coefficient uncertainty bounds established before probe entry are also 
shown in the figure. Note that the even in the extreme case of assuming all errors are in 
the aerodynamics (the LN case) the axial force coefficient is within the pre-entry 
boundaries. As expected, the LN case has more variation in the axial coefficient than the 
HN case. Also, the HN results (which is closer to the nominal aerodynamics) are roughly 
the mean of the LN case as would also be expected since the HN case does not “chase” 
the measurement data points as rapidly as the LN case does. While a bit higher towards 
the end of the entry phase (around 1000 m/s), this reconstructed case shows little 
adjustment in Ca is required beyond the nominal predicted values to be consistent with 
the flight data. This result again indicates that the pre-entry aerodynamics predictions 
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appear to have been accurate within the degree of uncertainty associated with the 
analysis. 


Parachute Phase Assessment 

For the entry phase, the dominant aerodynamic force is acting in the direction of 
the one accelerometer (i.e., axially). The parachute phase was more difficult to 
reconstruct for several of reasons. First, while the assumption was that the parachute was 
only influenced by drag, some side force would be created when the vehicle velocity (or 
parachute) is not aligned with the vehicle axis of revolution (the measurement direction 
for the axial accelerometer). Without lateral accelerations, off nominal conditions (such 
as wind gusts) are difficult to verify. Also, there are multiple possible solutions that can 
match the axial acceleration, but not satisfy other measured constraints (such as total 
descent time or altitude rate of change) due to the unknown amount of lateral (and thus 
total) acceleration sensed by the system. Second, the probe spent over two hours on the 
drogue parachute at terminal velocity (very little acceleration). Finally, the acceleration 
data collected had nearly +/- 0. 1 m/s 2 of noise in the signal. 

Figure 17 shows the comparison of accelerations from a simulation only run, 
accelerometer bias reconstruction, unbiased flight accelerometer data, and gravity 
acceleration during the drogue parachute phases. Since the parachute should be at or very 
near terminal velocity, the measured acceleration should equal the acceleration due to 
gravity. Note that the simulation-only run and gravity have nearly the same acceleration 
after about 4000 sec; whereas, the unbiased accelerations are notably below the gravity 
values. A reconstruction case was run to determine the bias required to bring the 
measured accelerations inline with the gravity. Figure 18 shows the bias solution 
determined from the reconstruction run. The biased accelerations are shown in Fig. 17. It 
appears that the accelerometer is biased about 0.022 m/s 2 during this long drogue 
parachute phase. An accelerometer bias of this magnitude (2200 p-g) is unusually large 
and not anticipated for a space quality instrument, particularly since just prior to 
atmosphere entry the bias was determined to be about 23 p-g, as mentioned earlier. This 
change in bias by a factor of about 100 means that either the events occurring during the 
EDL, or the environment at Titan can substantially modify navigation sensors, something 
of which future planners need to be aware. 



Accel Bias (m/s ) Longitudinal Acceleration (m/s 2 ) 
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1 7 Accelerations and Accelerometer Data in Drogue Phase 



Figure 18 Estimated Bias to Measured Accelerometer Data 
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Also noted in Fig. 18 is a large variation in estimated bias prior to 4000 sec. This 
variation is also seen in the measured accelerations (see Fig. 5 between 2200 and 4000 
sec) and would appear to be some (as yet undetermined) phenomenon that occurred to the 
probe. One possible explanation would be a wind shear, but no solutions thus far have 
been consistent across all measured data. That is, a solution can be found which matches 
the axial acceleration profile using vertical wind shears, however this solution provides 
altitude rates inconsistent with those derived from the pressure and temperature 
measurements and a shorter duration descent than was recorded from probe radio 
transmissions and surface impact measurements. Thus, a bias prior to about 3800 sec 
was not established or used for solutions. 

State (Position, Velocity) Estimate Runs 

These trajectory reconstruction runs, focused on estimating the vehicle state 
(position and velocity) throughout the entire EDL trajectory, were done using the initial 
state determined from the Monte Carlo analyses as indicated above. The bias determined 
in the preceding section was also applied to these cases. The deceleration pulse during 
entry for this reconstructed trajectory is shown in Fig. 19. The region of maximum 
deceleration in this figure shows that while the reconstruction run captured the 
acceleration well, it did not overly adjust the state to meet the measurements. The region 
just after parachute deploy also shows good agreement between the flight and 
reconstructed accelerations. 



Figure 19. Huygens Flight Data Compared to State Estimate Reconstruction Run- Entry 
Phase 
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Figures 20 and 21 focus on the flight data comparison at the end of the Main 
parachute and throughout the Drogue parachute phases. These figures illustrate the 
agreement between the flight axial accelerations and those generated using the 
reconstructed trajectory. The reconstructed trajectory data matches well through the 
fifteen minute Main parachute phase and through most of the Drogue phase. The curves 
separate slightly between 2000 and 2500 sec, but come back into agreement as time 
progresses. Note that the reconstruction case follows the biased accelerometer data as 
desired (see Fig. 21). 



Figure 20. Comparison of Flight Data to State Estimate Reconstruction Run - Main and 
Drogue Parachute Phases 
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Figure 21. Comparison of Flight Data to State Estimate Reconstruction Run - Drogue 
Parachute Phase 

A comparison of the altitude versus time for the reconstruction case and the radar 
altimetry data for unit A is included in Fig 12. As noted in the figure, the reconstruction 
case is slightly below the radar altimetry, by about 5 km, at the highest altitude but then 
reduces to less than 100 m by touchdown. A comparison of several reconstruction cases 
is shown in Fig. 22. As noted in this figure, ah of the cases show a similar comparison to 
the radar data. This result would indicate that a scale factor should be applied to the radar 
data. 
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Figure 22. Altitude Profile for State Estimate of Several Cases -Radar Altimetry 
Comparison 

Figure 23 shows the effect of a 0.93 scale factor on the radar data when compared 
to the same cases as shown in Fig. 22. This very good agreement between adjusted flight 
data and multiple reconstructed trajectories suggests that an adjustment of the radar data 
may be required. Alternatively, other possibilities exist that would also affect the radar 
measurements, namely probe tilt and radar electronics temperature sensitivity, but those 
possibilities were not considered this analysis. 
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Figure 23. Altitude Profile for State Estimate of Several Cases-Radar Altimetry Scale 
Factor 


CONCLUSIONS 

Trajectory reconstruction of the Fluygens probe entry and descent at Titan was 
completed using two separate approaches: a traditional method and a POST2 simulation 
based method. Both reconstruction methods showed agreement with the DTWG and 
Carlo Bettanini generated altitude profiles that are in conflict with the altimetry flight 
data returned from both radars. A scale factor of about 0.93 would bring the radar 
measurements into compliance. Further analysis shows that the Titan-GRAM atmosphere 
model was within 1% of the reconstructed density below about 160 km altitude (or where 
the onboard pressure and temperature measurements began). A bias in the accelerometer 
measurements of around 0.022 m/s 2 after about 3800s from atmospheric interface was 
also suggested from the data analysis. Finally, an assessment of the entry aerodynamics 
model (generated prior to entry) showed that the variation in the aerodynamics as shown 
by the flight data was well within the 3 -sigma bounds established by the aerodynamicists 
before entry. 


ACKNOWLEDGEMENTS 

The authors would like to thank the following groups and individuals: NASA 
Engineering Safety Center for funding this effort; the Huygens project scientists for 


AAS 07-226 


sharing their data; Jean-Pierre Lebreton of ESA for his leadership of the Huygens mission 
and allowing our involvement; Carlo Bettanini of University of Padova, Italy, David 
Atkinson of the University of Idaho, and Bobby Kazeminejad of the German Space 
Operations Center (DLR) for their data and discussions regarding their reconstruction of 
the Huygens probe trajectory. 


REFERENCES 

1 Lebreton, J.-P., et al., “An overview of the descent and landing of the Huygens probe on 
Titan”, Nature, Vol 438, 8 Dec 2005, pp 758-764. 

2 Lebreton, J.-P., and D. L. Matson (1997), “The Huygens Probe: Science, payload and 
mission overview”, Space Science Reviews, Vol. 104, 2002, pp. 59-100. 

3 Kazeminejad, B., et al., “Huygens’ Entry and Descent through Titan’s Atmosphere - 
Methodology and Results of the Trajectory Reconstruction”, Planetary and Space 
Science, in Press, 2007. 

4 Bird, M. K. et al., “The vertical profile of winds on Titan”. Nature, Vol 438, 8 Dec 
2005, pp 800-802. 

3 Lulchignoni, M. et al., “In situ measurements of the physical characteristics of Titan’s 
environment”, Nature, Vol 438, 8 Dec 2005, pp 785-791. 

6 Tomasko, M. G. et al. Rain, winds and haze during the Huygens probe’s descent to 
Titan’s surface.” Nature, Vol 438, 8 Dec 2005, pp 765-778. 

8 Kazeminejad, B., Methodology > Development for the Reconstruction of the ESA Huygens 
Probe Entry and Descent Trajectory. Ph.D. thesis, Karl-Lranzens Univ., Graz, Austria, 
2004. 

9 Lolkner, W.M., et al., “Winds on Titan from ground-based tracking of the Huygens 
probe”, Journal of Geophysiccd Research, Vol. Ill, 2006,E07S02. 

10 NASA Engineering and Safety Center, “Assessment of the Cassini/Huygens Probe 
Entry, Descent and Landing (EDL) at Titan,” NESC-RP-05-67, May 2005. 

11 Striepe, S.A., et al., “Program to Optimize Simulated Trajectories (POST II), Vol. II 
Utilization Manual.” Version 1.1. 6.G, January 2004, NASA Langley Research Center, 
Hampton, VA. 

12 Spencer, D.A., Blanchard, R.C., Braun, R.D., Kallemeyn, P.H., and Thurman, S.W., 
“Mars Pathfinder Entry, Descent, and Landing Reconstruction,” Journal of Spacecraft 
and Rockets, Vol. 36, No. 3, 1999, pp. 357-366. 


27 



13 Tolson, R.H., et al., “Application of Accelerometer Data to Mars Odyssey Aerobraking 
and Atmospheric Modeling,” Journal of Spacecraft and Rockets , Vol. 42, No. 3, 2005, 
pp. 435-443. 

14 Witkowski, A., et al., “Mars Exploration Rover Parachute System Performance”, 

AIAA Paper 2005-1605, 18 th AIAA Aerodynamic Decelerator Systems Technology 
Conference and Seminar, Munich, Germany, May 2005. 

15 Spencer, D.A. , Braun, R.D., “Mars Pathfinder atmospheric entry - Trajectory design 
and dispersion analysis,” Journal of Spacecraft and Rockets (ISSN 0022-4650), 
Sept./Oct. 1996, vol. 33, no. 5, pp. 670-676. 

1(1 Desai, P.N., Schoenenberger, M., and Cheatwood, F.M., “Mars Exploration Rover Six- 
Degree-Of-Freedom Entry Trajectory Analysis,” Paper AAS 03-642, AAS/AIAA 
Astrodynamics Specialists Conference, Big Sky, Montana, August 3-7, 2003. 

17 Desai, P.N., Lyons, D.T., “Entry, Descent, and Landing Operations Analysis for the 
Genesis Re-Entry Capsule,” Paper AAS 05-121, 15 th AAS/AIAA Space Flight 
Mechanics Conference, Copper Mountain, CO, January 23-27, 2005. 

Ix Desai, P.N., Mitcheltree, R.A., and Cheatwood, F.M., “Entry Trajectory Issues for the 
Stardust Sample Return Capsule,” International Symposium on Atmospheric Reentry 
Vehicles and Systems, March 16-18, 1999, Arcachon, France. 

19 Striepe, S.A., Way, D.W., Dwyer, A.M., and Balaram, J., “Mars Smart Lander 
Simulations for Entry, Descent, and Landing,” AIAA Paper 2002-4412, AIAA 
Atmospheric Flight Mechanics Conference, Monterey, CA, August 2002. 

20 Justus, C.G., et al, “Connecting atmospheric science and atmospheric models for 
aerocapture at Titan and the outer planets,” Planetary and Space Science. Surfaces and 
Atmospheres of the Outer Planets, Their Satellites and Ring Systems (ISSN 0032-0633), 
April 2005, Volume 53, no. 5, pp. 601-605. 

21 Gelb, A. (ed.). Applied Optimal Estimation, MIT Press, Cambridge, MA, 1974. 

22 Brown, R.G, and Hwang, P.Y.C., Introduction to Random Signals and Applied Kalman 
Filtering, John Wiley and Sons, New York, NY, 1992. 



